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A  AS  VkA.'T 

An  undorsta  nd  i  rvi  o  •'  the  mechanical  oroporLios  anil  beha¬ 
vior  of  the:  intr-rvertobra  1  disc  is  critical  to  several 
current  areas  of  research .  A.nonq  these  are  the  study  of  the 
effects  of  extreme  gravitational  forces  on  the  air  crews  of 
h  iqh  oer  forma  noc*  aircraft,  the  related  problem  of  soinal 
injuries  due  to  aircraft  e  lection ,  and  the  study  of  disc 
leqoneration.  This  study  was  undertaken  in  order  to 
construct  a  realistic  anilytic-ai  model  of  the  intervertebral 
joint  using  the  finite  e lament  approach. 

Experiments  have  shown  that  healthy  intervertebral 
joints  exhibit  creep  when  subjected  to  axial  load.  A  pre¬ 
vious  study  of  this  time  dependent  behavior  employed  i 
simplified  visco-o  last ie  model  of  the  disc  which  neglected 
its  inhomopeneity .  Thouah  the  homoqe neous  mode  1  success¬ 
fully  simulated  the  ext  trnally  observed  response  of  its 
joint,  it  could  not  adequately  portray  the  internal  response 
of  the  disc:,  especially  the  interaction  of  the  annulus 
fibrosis  and  the  nucleus  pu l port in. 

The  focus  or  this  invest iqat ion  is  the  behavior  and 
interaction  of  the  nucleus  and  annulus  under  loadinq, 
notinq,  in  particular,  the  role  of  the  annulus  as  a 
restrainin'-:  mechanism  on  the  nucleus.  An  ax  i  symmetric 
finite  element  mode l  is  employed  which  incorporates  a  linear 
visco-elastic  constitutive  relation  for  the  annulus  and 


nucleus.  This  relation  is  based  upon  a  three -parameter 
Kelvin  solid.  The  visco-elastic  constants  are  found  by 
matching  one  dimensional  axial  experimental  data  with  this 
two  dimensional  model.  The  nucleus,  which  is  composed  of  a 
series  of  concentric  lamellae  of  collagen  fibers,  has  been 
modelled  as  a  \7.  ply  struct  ire  of  orthotropic  material.  The 
cortical  bone,  traoecu  lar  bone,  and  th  >  bony  end  plate  are 
assumed  to  be  isotropic  in  i  l i noa  fly  elastic. 

Results  are  presented  which  depict  the  displacement 
profiles  and  stress  re  1  Lstr ibut Lons  occurring  as  a  con¬ 
sequence  of  the  interaction  if  the  annulus ,  nucleus,  and 
bony  enl  plate  under  axial  comorer.sive  loading.  Horizontal 
profiles  of  each  of  the  stress . components  are  presented  for 
three  regions  of  the  point;  the  trabecular  bone,  its  bony 
end  plate,  and  the  intervertebral  disc.  The  latter  pro £i Los 
clearly  show  the  importance  of  the  nucleus  as  a  load 
carrying  structure  and  the  action  of  the  annulus  in 
restraining  the  outward  flow  of  the  nucleus. 

Also  presented  are  the  variations  in  stress  with  time 
in  three  of  the  lamina  o c  the  annulus  (an  inner,  middle,  and 
outer  ply).  Results  indicate  that  the  orthotropic  material 
properties  of  the  annul  is  have  a  significant  impact  upon  its 
time  dependent  behavior  and  should  be  included  when 
modelling  the  point  . 
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A  FINITH  KbF.Mf’NT  ANALYSIS  OF  THF 
V I  SCO- R I. ART  l  C  [NTMRA'TION  OF  THR  NUCbRNS 
PULPOSIS  AMO  ANMUT.UR  FIBROSIS  WITHIN  THF 
HUMAN  IMTKRVHRTRBRAL  JOINT 

I .  INTRODUCTION 

1 . I  Purpose 

The  purpose  of  this  thesis  is  to  construct  a  realistic 
analytical  mold  of  the  human  intervertebral  joint,  using  the 
finite  element  approach.  It  is  an  extension  of  an  earlier 
work  (Ref  10)  which  use-?  an  a  x  isymmot  r  ic ,  visco-elastic 
model  of  the  disc  to  study  its  dm*'  depots  lent  behavior. 

This  the si s  further  refines  that  model  by  accounting  for  the 
actual  physical  structure  of  the  annulus  fibrosis  as  well  as 
its  orthotropic  properties.  ~’ho  focus  of  this  investigation 
Is  the  behavior  rind  inter  iction  of  the  nucleus  pulposis  and 
annulus  fibrosis  under  axial  compressive  loading.  The 
stress  redistributions  and  displacement  fields  occurring  as 
a  result  of  the  inclusion  of  a  separate,  visco-elastic  annu¬ 
lus  are  i  nves 1 1  ga t.ed .  Of  part  icular  interest  is  the  role  of 
the  nucleus  as  a  load  cirrying  structure  and  the  action  of 
the  annulus  as  as  restraining  mechanism  on  the  outward  flow 
of  the  nucleus  under  loading. 
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I  .  2  Background 

An  understanding  of  t.h->  mechanical  properties  and  be¬ 
havior  under  loading  of  the  human  intervertebral  joint  is  per¬ 
tinent  to  several  current  ureas  of  research.  Among  these 
ire  the  study  of  th->  effects  of  extreme  gravitational  forces 
on  the  aircrews  of  hinh  performance  aircraft,  the  related 
problem  of  spinal  injuries  following  aircraft  ejection,  the 
treatment  of  conqerital  spina!  deformities,  and  the  develop¬ 
ment  of  suitable  materials  for  disc  replacement. 

Although  a  sizeable  amount  of  data  has  been  gathered 
from  experimental,  testing  of  vertebral  segments,  the 
complexity  of  the  structure  has  prevented  the  accurate 
determination  of  material  property  constants  for  the  various 
regions  of  the  joint.  The  most  complicated  problems  in  anal¬ 
ysis  are  posed  by  the  annulus  fibrosis.  Composed  of  a 
series  of  concentric  lamella*'  of  collagen  fibers,  the 
orthotropic  nature  of  the  annulus  requires  that  at  least 
seven  distinct  material  constants  be  determined  for  a  linear 
analysis  (Ref  1).  Information  of  this  type  is  very  limited 
at  the  present  time. 

Experimental  investigations  into  the  load  deflection 
behavior  of  the  intervertebral  joint  were  carried  out  by 
Drown  (1957),  Nach^mson  (I960),  Mirsch  (1965),  and 
Rolander  (1966).  Galante  (1967)  examined  the  effects  of 
inhomogeneity  and  anisotropy  on  the  tensile  properties  of 
the  annulus  fibrosis  (Ref  7).  Markolf  (1972)  conducted 
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axial  compression,  and  tension  tests  on  thoraco Icmba r  discs 
and  noted  significant  differences  in  response  between  the 
lumbar  and  thoracic  .Uses  (Ref  17).  Kazarian,  et  al.  (  1 9 7 S ) 
reported  the  creep  characteristics  of  intervertebral  joints 
under  constant  axial  load  (Ref  11).  More  recently,  Kazarian 
and  Kaleps  (1979)  presented  a  method  for  determining  the 
■fount's  Modulus  and  viscosity  coefficient  for  a  thr>m- 
parameter-sol  id  model  based  on  their  earlier  data  (Ref  14). 
Rums  (1980)  developed  a  mathemat i  ca  1  analysis  technique  for 
determining  unique  parameter  values  Cor  the  three-parameter- 
solid  model  using  experimental  strain  data. 

Finite  element  analysis  techniques  have  !>een  success- 
rul!y  employed  by  several  researchers  to  model  the  joint. 
Mpilker,  Be ly tschko,  and  Kulak  used  axisynmotric  finite  ele¬ 
ment  models  to  study  tin.*  time  independent  behavior  of  the 
disc  (Refs  21,  1,  15).  Re ly tschko 1  s ,  et  al.  model  (1974) 
assumed  an  incompressible  nucleus  in  a  hydrostatic  state  of 
stress  with  an  orthotropic,  inhomogeneous,  annulus  fibrosis. 
This  study  led  to  the  conclusion  that  the  material  ani¬ 
sotropy  of  the  disc  could  not  be  omitted  from  the  analysis 
without  causing  large  errors.  It  was  also  found  that  a 
linear  analysis  could  not  account  for  the  differences  be¬ 
tween  thoracic  and  lumbar  disc  behavior  or  the  differences 
between  tensile  nnd  compressive  response  (Ref  15).  Kulak, 
et  al  (1976)  expanded  on  the  earlier  work,  using  a 
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nonlinear,  elastic  constitutive  relation  for  the  annulus  to 
study  the  aforementioned  differences. 

Most  recently,  Hinrichson  (1980)  employed  an  axisym- 
motric  finite  element  model  which  included  material  inhomo- 
genoity  and  visco-elasticity  to  study  the  creep  response 
reported  by  Kazarian  (Ref  10).  The  disc  was  idealized  as  a 
single,  homogeneous,  isotropic,  linear,  visco-elastic 
substance.  This  model  successfully  reproduced  the  exter¬ 
nally  observed  response  of  the  joint  but  due  to  its 
simplifications,  could  not  adequately  portray  the  internal 
respons'  of  the  disc,  in  particular,  the  interaction  of  the 
annulus  fibrosis  and  the  nucleus  oulposis. 

1.3  Anatomy  of  the  Intervertebral  Joijut 

The  human  spine  is  a  flexible  column  composed  of  a 
series  of  24  segmented  bones  called  the  vertebrae  (see 
Fig  1.3-A).  Fach  vertebra  consists  o':  a  ventral  body  and  a 
set  of  posterior  elements  (Fig  1.3-3).  One  of  those  ele¬ 
ments  is  a  spinous  projection,  or  process,  which  extends 
downward  and  backward.  These  processes  may  be  felt  as  a 
series  of  bumps  down  the  back.  In  addition,  there  are  two 
transverse  processes,  one  to  either  side,  which  provide 
attachment  for  muscles  and  ligaments.  The  vertebral  body 
itself  is  composed  of  soft,  trabecular  bone  encased  in  a 
thin  shell  of  relatively  hard  cortical  bone.  On  the  upper 
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Human  Vertebral  Column. 


( from  Cunningham's  Textbook  of  Anatomy  by  G.  J.  Romanes, 
Univorsitv  of  Ox  fort''  Press) 
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and  lower  surfaces  of  the  body  are  thin  plates  of  cortical 
bone  called  bony  end  plates. 

Between  each  pair  of  vertebrae  is  a  cartilaginous 
i ntervertebra L  disc.  The  discs  vary  in  size,  shape,  and 
thickness  at  different  spinal  levels.  The  lumbar  discs  are 
the  thickest,  while  those  in  the  thoracic  region  are  rela¬ 
tively  thin.  The  disc  Is  composed  of  a  soft  center  called 
the  nucleus  pulposis  which  is  encircled  by  a  tough, 
flexible,  f ibrocarta l ig i nous  ring  called  the  annular 
f i Pros i s . 

The  nucleus  pulposis  is  composed  of  a  collagenous  lat¬ 
tice  enmeshed  in  s  mueoorotein  gel.  Its  water  content  ranges 
from  88  percent  at  birth  to  approximately  69  percent  in 
middle  age.  experiments  carried  out  by  Nachemson  (Ref  18) 
indicate  that  the  nucleus  in  normal  discs  behaves 
hydrostatical  ly . 

The  annulus  fibrosis  is  formed  by  a  series  of  con¬ 
centric  encircling  lamellae  (see  Fig  1.3-C).  The  .lamellae 
consist  of  collagen  fibers  which  have  two  well  defined  areas 
of  orientation.  The  fibers  run  in  a  single  direction, 
alternating  from  the  previous  one  and  aligned  at  approxi¬ 
mately  30°  to  the  horizontal  (Ref  1). 

Fxterior  to  the  disc  proper  are  two  ligaments?  the 
anterior  and  posterior  longitudinal  ligaments.  Although 
these  ligaments  transmit  considerable  force  during  complex 
loading  (Ref  16),  they  do  not  appear  to  be  significant  as 
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load  carrying  structures  under  axial  loading  only. 

Therefore,  in  order  to  focus  on  the  disc  interaction,  the 
ligaments  are  usually  removed  during  axi3l  testing  of  the 
disc. 

1.4  Genera  1  Approach  a_nd  Assumptions 

The  vertebra  is  modelled  as  an  axisymmetric  body  which 
is  symmetric  with  respect  to  th->  centerline  of  the  interver¬ 
tebral  disc.  Only  the  body  itself  is  modelled,  based  upon 
the  assumption  that  little  if  »ny  Load  is  transmitted  by  the 
anterior  and  posterior  longitudinal  ligaments  during  axial 
loading  (Ref  l).  Pictured  in  Fig  1.4-A  are  two  adjacent  ver¬ 
tebrae  and  their  intervertebral  disc.  The  dotted  lines 
indicate  the  truncated  portion  of  the  joint  used  in  the 
finite  element  model.  Fig  I. 4-B  is  an  enlarged  diagram  of 
the  model  showing  the  various  regions  of  the  joint.  The 
outer  radius  of  the  nucleus  pulposis  is  taken  to  be  0.707 
times  the  outer  radius  of  the  disc  (Ref  1). 

Kazarian's  research  (Ref  Id)  indicates  that  the  ver¬ 
tebral  body  alone  does  not  exhibit  creep  under  constant 
axial  load.  Therefore,  it  is  assumed  that  the  trabecular 
bone  in  the  vertebral  core  and  the  cortical  bone  in  the 
outer  shell  and  end  plates  are  linear,  elastic  materials. 

Since  the  vertebral  body  does  not  appear  to  contribute 
to  the  observed  creep  behavior  of  the  joint,  it  is  assumed 
that  the  disc  is  the  visco-elastic  medium.  However,  at  the 
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time  of  this  analysis,  it  was  unknown  whether  the  visco¬ 
elastic  response  should  be  attributed  to  the  nucleus,  the 
annulus,  or  both.  As  a  first  quess,  it  was  assumed  that  the 
nucleus  alone  was  the  visco-elastic  medium.  The  assumption 
was  based  upon  Kazarian’s  {Ref  14)  work  with  degenerated 
discs.  The  nuclei  of  these  discs  exhibited  varying  degrees 
of  herniation  and/or  desiccation,  mainly  due  to  ageing. 
Kazarian  reported  that,  unlike  healthy  discs,  the  degenerated 
discs  did  not  exhibit  creep  response.  Therefore,  the  annu¬ 
lus  was  originally  assumed  to  be  linearly  elastic.  During 
the  course  of  the  investigation,  however,  it  was  decided  to 
model  both  the  nucleus  and  the  annulus  visco-e 1 ast ica lly . 

The  steps  which  led  to  this  decision  are  discussed  in  detail 
in  Chapter  IV. 

An  existing  finite  element  program  was  adopted  for  use 
in  this  investigation.  The  original  program  was  written  by 
Hinnerichs  (Ref  0)  to  handle  plane-stress,  plane-strain 
problems  in  isotropic  materials.  Hinrichsen  (Ref  10) 
modified  the  program  to  accommodate  axisymmetric  structures 
and  to  account  for  visco-elastic  response  in  inhomogeneous 
materials.  In  order  to  account  for  the  orthotropic  proper¬ 
ties  of  the  annulus,  the  program  was  again  modified. 

In  the  following  chapter,  the  basic  theory  behind  the 
finite  element  analysis  method  is  reviewed.  Included  there 
is  the  derivation  of  the  orthotropic  material  property 
matrix  used  to  modify  the  existing  program.  The  procedures 
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used  to  validate  the  modified  program  are  described  in 
Chapter  III.  Chapter  IV  Inscribes  the  details  of  the 
investigation,  including  selection  of  mesh  size  and  deter¬ 
mination  of  material  constants.  The  final  results  are  pre¬ 
sented  in  the  form  of  displacement  profiles  and  stress 
re distri buttons. 
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THEORY 


r  i . 

In  the  followin'.}  sections,  the  basic  concepts  and 
equations  used  in  finite  element  analysis  are  briefly 
reviewed.  The  derivation  of  the  orthotropic  material  prop¬ 
erty  matrix  used  to  modify  the  finite  element  program  is 
discussed  in  Section  2.3.  In  the  final  section,  an  overview 
of  linear  visco-elastic  theory  and  the  procedure  by  which  it 
is  incorporated  into  finite  element  analysis  is  presented. 

2.1  F inite  Element  Method 

The  basic  concept  of  the  finite  element  method  when 
applied  to  problems  of  structural  analysis  is  that  a  complex 
structure  can  be  subdivided  into  a  finite  number  of  discrete 
elements,  in  each  of  which  the  stresses  or  displacements  can 
be  represented  by  relatively  simple  functions.  The  number 
of  elements  in  the  structure  must  be  large  enough  so  that 
the  displacement  function  for  each  element  closely  approxi¬ 
mates  the  exact  displacements  in  that  particular  region. 

The  solution  of  the  complete  system  is  then  represented  by 
an  assembly  of  its  elements.  As  the  size  of  the  elements 
becomes  small,  the  solution  should  converge  to  the  exact 
solution  for  the  structure. 

The  displacement  function  which  is  chosen  to  represent 
the  behavior  of  the  elements  must  meet  three  conditions  to 
assure  solution  con vergence .  First,  displacements  must  be 
continuous  l>etween  adjacent  elements.  Second,  rigid  body 


1  ispi  acomonts  uuy  o  >t  cause  st.raininq  of  the  element. 
Finally,  as  the'  element  size  approaches  zero,  Lae  strain 
must  reach  a  constant  value.  The  simplest  displacement 
function  which  satis! ios  these  conditions  is  a  linear 
polynomial.  A  polynomial  of  infinite  order  corresponds  to 
the  exact  solution.  However,  as  the  element  size  becomes 
sra!  I  ,  a  finite  polynomial  provides  a  simple  approximation 
to  the  exact  solution. 

Figure  2.  1-A  snows  a  typical  trianqular  element  with 
nodes  1,  2,  3  numbered  in  »  counterclockwise  direct  ion. 


Pic.  2.1 -A. 

The  behavior  of  t.h.o  element  is  described  by  six  deqrees  of 
freedom  which  cor  respond  to  the  nodal  displacements: 


{  u]  vi  it  2  v  2  u3  Vj  )  • 


(2-1) 


The  linear  polynomial  which  defines  the  displacements  within 
the  element  is 


(2-2) 
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This  may  also  tie  written  in  matrix  form  as: 
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Fly  inversion  and  Kick-substitution  into  Eq  2-2  the  displace¬ 
ment  field  may  now  be  given  in  terms  of  the  nodal  displace¬ 
ments  as 
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and  A  is  the  area  of  the  t r i a n g u 1 a r  element  given  by: 
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2.2  Stress-Strain  Eolations  for  an  Ax^isymme tr ic  Bodjy 

The  total  strains  at  my  yioint  within  an  element  are 
defined  in  terms  of  tfv  d i so  1 acement  derivatives  as 
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by  appropt  i  a  t  o  Iv  <1:  f  f  ot'-nt .  i  a  t  inq  Eqs  2-4(a)  and  2-4fb), 
the  st  rain  -.i  i  sp  la  cr*  r.o  p  t  aquation  (7)  may  ho  rewritten  in  the 
following  form 
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with  respect  to  r. 

etc. 

For  a 

linear , 

elastic,  axisymmetric 

body,  the  stress- 

strain  relationship 

may 

he  expressed  by 
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(2-10) 
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where  {>,0}  represents  any  initial  strains  and  [D]  is  the 
material  property  matrix. 
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2.3  Derivation  of  the  Or  t  hot  r  op  ic  Material  Property  Matrix 
For  isotropic  ma  t>->r  i  a  Is  such  as  the  cortical  and  trabe¬ 
cular  hone,  nucleus  ;>u !  nos  i  s ,  and  Irony  end  plate,  the 
material  property  matrix  takes  the  following  relatively 
simple  form: 
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(2-11) 


where  v  is  Poisson's  ratio  and  R  is  the  bulk  modulus  of  the 
material . 

For  orthotropic  materials  whose  principal  directions 
coincide  with  the  axis  system  the  material  property  matrix 


becomes 
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G  i  2 


shear  modulus  in  the  1-2  plane 


vij  =  Pois. son's  ratio  for  transverse  strain  in  the 
j-dir;ection  when  stressed  in  the  i-direction. 

The  annulus  fibrosis,  however,  is  an  orthotropic 
material  whose  principal  directions  do  not  coincide  with  the 
axisymmetric  body  coordinate  of  the  model.  It  is  composed 
of  a  series  of  concentric  lamellae  in  each  of  which  the 
fibers  run  in  a  single  direction,  alternating  from  the  pre¬ 
vious  one  and  aligned  at  an  angle  $  to  the  0  axis.  Each 
lamella  is  orthotropic  in  its  plane  with  respect  to  coor¬ 
dinates  normal  and  tangent  to  the  fiber  direction. 

Pig  2.3-A  shows  the  relationship  of  the  principal  directions 
to  the  body  coordinates  in  a  segment  of  a  single  lamella. 


Fig  2.3-A.  Coordinate  Systems. 

The  stresses  and  strains  in  the  principal  directions 
are  related  to  those  in  the  axisymmetric  body  coordinates  by 
the  transformation  matrix  [ T )  as  follows. 
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Therefore,  combining  Cqs  2-13, 
noting  the  general  relationship  for 
the  principal  axis  direction, 
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2-14,  and  2-15,  and 
stresses  and  strains  in 
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we  can  now  derive  the  ( D ’ j  matrix  which  relates  stresses  and 
strains  in  the  body  axes  as  follows: 


"  Of)  " 

-  03  - 

0  2 

*  [T] 

a2 

)  i 

°r 

°3 

t  °rz 

^°12  . 

IT]  [D] 


'el 

I  G2  , 

E3 

Y12  _ 


18 


J 


ki — a-.,  u. 


[T]  ID]  [TJT 


E0 


or , 


°0 

e  0 

0 Z 

[O'  ! 

Ez 

> 

1  r  /. 

_Erz 

(2-17) 


where  [O'  )  =  [T]  (D)  [T]T 


(2-18) 


Note  that  in  assuming  that  r-axis  and  the  3-axis  am 


coincident,  we  have  neglected  the  physical  curvature  of  the 
annulus.  In  reality,  the  1-2  plane  o£  the  lamellae  is 


tilted  with  respect  to  the  0 -z  plane  by  a  small  angle  which 
varies  according  to  position  within  the  annulus.  Pig  2.  3-B 
is  an  exaggerated  diagram  of  the  annular  curvature  which 
points  out  the  angle  in  question.  In  order  to  account  for 
this  anqle  in  the  ( I) I  matrix,  a  second  transformation  would 
be  required.  However,  since  the  angle  is  small  (=  12°),  the 
cos2  and  sin2  terms  approach  1  and  0  respectively.  The 
transformation  matrix  ( Pq  2-18)  essentially  becomes  an  iden¬ 
tity  matrix.  Therefore,  for  the  purpose  of  this  analysis, 
the  angle  of  the  tilt  may  be  safely  ignored. 

Based  upon  the  foregoing  theory,  a  computer  subroutine 
which  calculates  the  [n'i  matrix  was  written  and  incor¬ 
porated  into  the  existing  axisymmetric  finite  element 
program.  The  subroutine,  and  the  methods  used  to  validate 
it,  are  described  in  detail  in  Chapter  III. 
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2.4  St i f f ness  Matrix 

The  general  equation  for  the  elastic  stiffness  of  an 
element  is 


IK  1  =  /  [  B 1 T  [  D  ]  [  B 1  <3  (volume)  (2-19) 

Vo  1.  ume 


Tor  an  axisymmetric  body,  the  volume  integral  must  he  taken 
over  the  entire  ring  of  material.  Summing  the  contributions 
from  each  of  the  elements,  the  stiffness  matrix  for  the 
structure  is  then 

N 

f  K !  =  l  2  n  /  [B]T  (D)  [B]  r  dr  dz  (2-20) 

n  =  1 

where  N  is  the  number  of  elomen-ts  in  the  structure. 

Since  the  [B]  matrix  contains  terms  which  depend  on  the 
coordinates,  straightforward  integration  will  not  suffice. 
However,  the  integration  may  be  approximated  by  evaluating 
[B]  at  a  centroidal  point.  The  centroid  is  given  by 


r  =  (rl  +  r2  +  r3)  /  3  (a) 

z  =  (ZL  +  z2  +  z  3 )  /  3  ( b ) 


(2-21) 


and  the  ( B ]  matrix  evaluated  at  the  centroid  is  denoted  by 
[Bj.  The  final  expression  for  the  elastic  stiffness  matrix 
is  then 


N 

[K]  =2  2  n  A  [BIT  [D]  [B]  r  (2-22) 

n -  .1 


with  [D]  given  by  Eqs  2-18,  2-12,  or  2-11,  depending  upon  the 
material  and  its  orientation. 
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2.5  Noda 1  Force  Vector 

The  external  nodal  forces  represent  a  combined  effect 
of  the  force  acting  along  the  whole  circumference  of  the 
circle  forming  the  element  "node”  (Ref  24).  For  the  pur¬ 
poses  of  this  analysis,  these  forces  shall  simply  be  denoted 
by  the  vector  {P}. 

The  nodal  forces  duo  to  initial  strain,  which  are  later 
used  in  the  solution  of  the  visco-elastic  problem,  may  be 
expressed  as 

N 

'Fef  =  >:  /  [B)'r  [D]  {e0}  r  dr  dz  (2-2  3) 

n=l 

The  same  approximation  procedure  which  was  used  to 
integrate  the  stiffness  equation  may  be  employed  here.  The 
final  result  is 

N 

{ F° }  -  2  2 w  A  [  B 1  [ 0 ]  { e0 }  r  (2-24) 

n  - 1 

2 . 6  Elastic  Solution 

Having  assembled  the  stiffness  matrix  and  nodal  force 
vectors,  we  may  solve  for  the  displacements  using  the  prin¬ 
ciple  of  virtual  work.  The  external  work  done  by  the  nodal 
forces  is  equated  to  the  internal  work  done  by  the  stresses 

to  yield  the  general  relation 

[K]  tu}  =  (P)  +  (FS)  (2-25) 

This  is  the  governing  equation  for  the  elastic  response 
of  a  discretized  structure. 


22 


2 . 7  Visco-El a s t i c  So 1 u t ion 

lip  to  this  point,  wo  have  dealt  only  with  elastic 

materials  Cor  which,  according  to  Hooke's  Law,  stress  is 

always  directly  proportional  to  strain  but  independent  of 

strain  rate.  Visco-elastic  materials,  like  the  nucleus 

pulposis,  behave  very  differently.  These  materials  display 

creep,  that  is,  an  increasing  deformation  under  sustained 

load,  with  the  rate  of  strain  depending  on  the  stress. 

The  visco-elastic  behavior  of  materials  may  be  modelled 

by  various  combinations  of  two  .simple  elements;  the  spring 

and  dashpot .  The  constitutive  equation  for  a  helical  spring 

which  obeys  Hooke's  Law  is 

o  =  ,E  e  (2-26) 

Similarly  Cor  the  dashpot 

o  =  n  e  (2-27) 

where  n  is  the  dashpot  constant  and  (°)  denotes  time 
d  i.  f  ferent  iat  ion  . 

The  simplest  combination  of  these  elements  is  the 
Maxwell  fluid.  The  elements  are  connected  in  series  and  the 
total  elongation,  c ,  is  a  sum  of  the  elongations  in  the 
spring  and  dashpot  as  shown  in  Fig  2.7-A. 

*  c  g 


Fig  2.7-A. 

2  3 


The  constitutive  relation  for  the  Maxwell  fluid  is 


(2-28  ) 


If  the  spring  and  dashpot  are  connected  in  parallel, 
the  resulting  model  is  called  a  Kelvin  solid: 


1 - VV'v- 

C) 


Fig  2.7-B. 


In  this  case,  the  total  force  a  is  a  sum  of  the  spring 
and  dashpot  forces.  The  resulting  constitutive  relation  for 
the  Kelvin  solid  is 


a  ~  Ft  +  ne 


(2-29) 


The  visco-elastic  model  used  in  this  analysis  is  a 
three-parameter  solid  composed  of  an  elastic  spring  coupled 
with  a  Kelvin  solid,  pictured  in  Fig  2.7-C. 


E 

- — 'W' - c 


Fig  2.7 -C . 


Since  the  instant  incous  elastic  deformation  is  depen¬ 
dent  only  on  the  spring,  it  can  be  separated  from  the  creep 
deformation  and  solved  by  the  finite  element  method.  If  the 
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constitutive  relation  for  the  three-parameter  solid  is 
integrated  over  time,  the  resulting  equation  is 


{  Af.c} 


[A]  { o  J  -  qg  ( rc ) 


(2-30) 


where  Aec  is  the  change  in  strain  due  to  creep  over  a  small 
time  interval  and  [A }  is  given  by 
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(  2-31) 


The  quantity  (Arc)  may  bo  considered  an  initial  imposed 
strain  {eg}.  Since  the  initial  strains  enter  into  the 
governing  f',q  2-2S  only  in  the  force  term  { Fe }  ,  a  method  of 
solution  to  the  visco-elastic  problem  now  presents  itself. 

As  time  is  increased  by  small  increments,  At,  the  governing 
equation  is  solved  at  each  time  step  for  the  values  of  {u}, 

{ a } ,  and  {(-c}.  The  values  for  (a)  and  {  ec  >  are  then  substi¬ 
tuted  into  Eq  2-30,  time  is  incremented,  and  the  process 
continues  until  a  maximum  time  is  reached. 

Further  discussion  of  the  application  of  numerical 
methods  to  the  solution  of  visco-elastic  problems  can  be 
found  in  publications  by  Minriehson  (Ref  10)  and  Zienkiewicz 


(Ref  23). 
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1 1 1 .  PROGRAM  VALIDATION 

Several  checks  were  made  to  assure  that  the  modified 
[D]  matrix  subroutine  produced  the  correct  coefficients  for 
each  type  element.  In  the  following  section,  the  internal 
workings  of  the  subroutine  are  briefly  described.  The 
remaining  sections  discuss  the  various  methods  used  to  vali¬ 
date  the  program. 

3 . I  Subroutine  DMATRI X 

The  modified  subroutine  works  in  a  very  simple  manner. 
The  main  program  was  altered  so  that  when  the  data  for  each 
element  is  read  in,  a  "tag"  denoting  whether  the  element  is 
isotropic  or  orthotropic  is  read  in  as  well.  The  tag, 
called  ISO  in  the  program,  is  equal  to  zero  for  isotropic 
elements  and  one  for  orthotropic  elements.  In  the  DMATRIX 
subroutine,  where  the  material  property  matrix  .is  computed 
for  each  element,  the  value  of  ISO  determines  whether  an 
isotropic  or  orthotropic  [D]  is  computed  for  the  element. 

The  procedure  is  diagrammed  in  Fig  3.1A. 
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Input  Element  Properties 
E  i  ,  IS  cl 

ISO  >  0  __  Yes 

No 

Compute  ISO  (D] 

From  Eg  2-11 


Fig  3.1-A 

An  initial  check  was  made  on  the  subroutine  as  a 
separate  entity  to  ascertain  that  data  was  correctly  read  in 
and  that  the  matrix  multiplications  produced  the  correct 
results.  Various  values  for  the  material  constants,  anti 
vjj,  were  input  and  the  resulting  [D]  matrices  compared  with 
hand  calculations.  When  this  check  was  satisfied,  the 
subroutine  was  added  to  the  main  program. 


N/ 

Compute  Ortho  [D] 
From  Eq  2-18 
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3.2  Comparison  With  a  Known  Structure 

Using  a  riosh  nearly  identical  to  Hinrichson's  (Ref  10) 
final  mash,  a  run  was  made  with  the  original  program.  In 
effect,  this  reproduce  1  Hinrichson's  final  results;  i.e.,  a 
visco-elastic,  isotropic,  homogeneous  disc.  The  same  data 
was  then  run  in  the  modified  program  with  ISO  =  0  for  all 
elements  in  the  disc.  The  results  of  the  two  runs  were 
identical,  proving  that  the  new  program  had  cortectly  com¬ 
puted  the  isotropic  [D]  matrix. 

3.3  Internal  Check 

To  assure  that  the  ISO  tag  functions  properly  in  the 
program,  two  more  runs  were  made  with  Hinrichson's  mesh. 
Noting  that  if  the  angle  ?  is  equal  to  zero  and  the 
orthotropic  material  constants  ate  given  by 
c  =  :q  -  f2  -  >--3 

v  "  v  1 2  =  v  2  1  ■=  v  1  3  =  v  3  l  =  v  2  3  =  v  3  2 
then  the  [n]  matrix  computed  from  Fq  2-18  is  the  same  as 
that  from  Fq  2-11.  Therefore,  if  elements  in  a  given  region 
are  tagged  ISO  =1,  i  -  0,  the  results  should  be  identical  to 
those  for  which  Ur*  same  elements  are  tagged  ISO  =  0, 
regardless  of  the  value  of  the  angle  (given,  of  course, 
that  the  material  constants  are  the  same). 

For  the  first  run,  ISO  was  set  equal  to  zero  in  all  of 
the  elements  of  the  disc  region.  F^r  the  second  run,  the 
same  material  constants  were  used,  but  in  the  annular 
region,  ISO  was  set  ■  •quel  to  one  and  the  angle  4>  set  equal 


to  30°.  As  desire*],  the.se  two  runs  produced  identical 
results . 

3.4  Force  Kfjui  1  i br ium 

As  a  .line  1  check,  on  the  validity  of  the  program,  we 
make  use  o!  the  integral  relation  between  the  applied  force 
and  the  resultant  stress, 

Pz  =  (  oz  uA  (3-1) 

A  ^ 

This  relation  must  hold  for  any  horizontal  i:ross  section  of 
the  body.  Fig  3. 4-A  shows  a  side  view  of  the  ioint  model 
through  which  a  horizontal,  slice  has  been  cut. 


Fin  3. 4-A.  Side  View  of  Loaded  Joint. 


If  the  in  stresses  along  the  horizontal  slice  are 
plotted  versus  the  radius,  then  the  area  under  the  curve, 
adjusted  by  a  factor  of  2nr,  must  equal  the  applied  force, 

P. 
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Applying  this  test  to  the  results  of  the  final  model  of 
the  disc,  a  stress  profile  is  taken  horizontally  through  the 
center  of  the  disc  at  t  -  0  seconds.  The:  resulting 
o ■’  versus  r  curve  is  shown  in  Fig  3.4-B,  The  area  under  the 
curve  is  calculated  in  three  intervals;  A3 ,  A2 ,  and  A3. 
Multiplying  these  areas  by  2n  times  their  respective 
centroidal  radii,  and  summing,  we  approximate  the  value  of 
the  integral  in  Fq  3-1.  These  computations  are  shown  in 
Fig  3.4-B. 

Figure  3.4-B  shows  that  the  area  approximation  of  the 
stress  integral  (17.7  Kiloponds)  is  very  close  to  the  actual 
applied  load  of  18.0  Kiloponds.  Since  the  final  model 
included  orthotropic  properties'  in  the  annulus  and  isotropic 
properties  in  the  nucleus,  we  may  be  satisfied  that  the 
program  correctly  computes  both  types  of  [Dj  matrices. 
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I v .  APPL I_C AT  ION  TO  THE  INTERVERTEBRAL  JOINT 

Kazarian's  (Ref  14)  experiments  involved  the  testing  of 
several  types  of  intervertebral  joint  specimens.  These  spec¬ 
imens  were  made  by  cutting  the  vertebral  bodies  parallel  to 
the  lines  of  demarcation  between  the  discs  and  the  bony  end 
plate  at  levels  of  maximum  vertebral  waisting.  The  speci¬ 
mens  were  then  measured  and  placed  in  the  test  apparatus. 

Each  was  subjected  to  a  dead  weight  load  of  1R.0  Kp 
(Kiloponds)  for  a  period  of  time  and  the  deflection  of  the 
top  edge  was  measured  (Ref  10). 

Since  the  objective  of  this  investigation  is  to  deter¬ 
mine  how  the  inclusion  of  orthotropic  material  properties  in 
the  annulus  affects  the  time  dependent  behavior  of  the  disc, 
wo  wanted  to  bo  able  to  compare  our  results  to  those  of 
Hinrichsen  (Ref  .1.0)  (i.e.,  to  an  isotropic,  homogeneous 

disc).  Therefore,  the  same  specimen  used  in  Hinrichsen's 
work  was  used  in  this  analysis;  specimen  number  65  (a 
T10-T11  segment).  So  that  comparisons  with  Kazarian's 
experimental  data  could  be  made,  the  same  1R  .  0  Kp  load  was 
applied  to  the  finite  element  model  as  was  applied  to  the 
test  specimen. 

4 . I  Determine tion  o f  M  a  ter i a  1  Constant s 

The  finite  element  model  is  divided  into  five  distinct 
regions  (refer  to  Fig  1.4-B).  The  vertebral  core  is  assumed 
to  be  an  isotropic,  homogeneous  material  with 
E  -  750  Kp/sq  cm  and  v  =  .25.  The  material  constants  for 
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the  cortical  bone  in  the  outer  shell  and  bony  end  plate  are 
taken  to  be  E  =  lf>0000  Kp/sq  cm  and  v  =  .25  (Refs  1,  10)  . 

The  annulus  is  orthotropic,  with  fibers  oriented  at 
+  t  and  -<f>  to  the  horizontal.  The  anqlo  <j>  is  assumed  to  be 
30°  (Refs  11,  l).  Assignment  of  the  orthotropic  material 

constants  is  based  upon  the  values  derived  by  Belytschko 
(Ref  1).  He  determined  the  values  of  E^ ,  E2,  and  Gi2  that 
best  matched  Galante's  (Ref  7)  experimental  results  at 
various  load  levels.  The  results  are  plotted  in  Fig  4.1-A. 
Taking  the  values  which  correspond  to  an  IB . 0  Kp  load,  the 
material  constants  are  as  follows:  E3  =  435  Kp/sq  cm, 

E2  =  10.9  Kp/Sq  cm,  and  G32  ~  7.2  Kp/sq  cm.  Noting  that 
each  lamella  is  orthotropic  in  its  plane  with  respect  to 
coordinates  normal  and  tangent  to  the  fiber  direction 
(Ref  1),  we  may  further  assume  that  F.g  is  equal  to  E2.  The 
Poisson  ratios  are  determined  from  the  relation 

vij  =  v j i  Ri  (4.1) 

Assuming  Pelytschko’s  value  of  0.45  for  v^,  we  have 

v  1 3  =  v12  =  ’45 

V31  ~  V21  =  *01128 


v23  =  v  32 


Assuming  that  the  material  is  compressible,  it  follows 

tha  t 

v23  lL  1  ~  v13  (4.2) 

with  the  eqjality  corresponding  to  incompressibility. 
Therefore,  the  value  of  V23  must  lie  less  than  0.55.  For  the 
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purpose  of  this  .analysis, 
assumed  to  be  0.'!'^. 


the  value  of  V23  (=  V32)  is 


In  the  nucleus,  a  range  of  values  for  the  Poisson's 
ratio  were  tried.  To  simulate  a  condition  close  to  an 
incompressible  nucleus,  the  value  of  v  was  set  equal  to 
0.<95.  An  intermediate  value  (the  same  value  assumed  by 
Hin-’ichsen )  was  set  at  0.48.  On  the  low  end,  v  was  set 
equal  to  .333  to  simulate  the  very  deformable  nature  of  the 
nucleus  material.  Noting  that  the  instantaneous  displace¬ 
ment  of  the  top  edqe  of  the  model  is  dependent  on  the  value 
of  the  Bulk  Modulus,  a  trial  and  error  method  is  used  to 
find  the  value  of  E  for  the  nucleus.  The  value  of  P.  is 
varied  untiL  the  displacement  of  the  model  at  t  ~  0  is  equal 
to  the  experimental  displacement .  It  was  found  Uiat  fhe 
experimental  displacement  of  0.014  cm  could  bo  matched  by 
the  model  using  any  value  of  v  with  an  appropriate  value  of 

R .  These  pairs  of  values  are  sho  below. 

v  =  .495  R  =  1.58  kp/sq  cm 

v  =  .48  E  =  6.2  kp/sq  cm 

v  -  .333  E  =  35.0  kp/sq  cm 

At  this  point,  only  the  visco-elastic  constants,  qp  and 
qi,  remain  to  be  determined.  The  procedure  is  similar  to 
the  method  described  above  for  finding  the  bulk  modulus.  In 
this  case,  the  displacement  versus  time  curve  of  the  top 
edge  of  the  model  is  compared  to  the  experimentally  deter- 
minod  curve.  Hinrichsen  found  that  the  slope  of  the  curve 
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is  controlled  by  the  value  of  qi,  while  the  qo  value  tends 
to  influence  the  maximum  displacement.  By  varying  the 
values  of  qo  and  qi  in  conjunction  with  one  another,  it  is 
possible  to  match  the  exact  shape  of  the  experimental  curve 
with  the  finite  element  model,  given  that  all  other  model 
assumptions  are  reasonably  accurate.  As  will  be  seen  in 
Section  4.4,  the  experimental  curve  cannot  be  matched  when 
the  annulus  is  assumed  to  be  elastic. 

4 . 2  Mesh  S i zes 

Four  different  mesh  sizes  were  employed  in  this 
analysis.  The  initial  mesh,  pictured  in  Fig  4.2-A,  is  very 
similar  to  the  mesh  developed  by  Hinrichscn  (Ref  10)  and 
satisfies  the  requirement  that  the  stresses  within  the  joint 
vary  smoothly  from  region  to  region.  The  mesh  contains  224 
elements  and  133  node  points.  Note  that  the  bottom  edge  of 
the  model  is  constrained  from  axial  displacement  and  that 
the  node  points  along  the  axis  of  rotational  symmetry  are 
constrained  from  radial  displacement.  As  discussed  in 
Chapter  III,  this  mesh  was  used  to  validate  the  modified 
proqram. 

A  coarser  mesh  (Fig  4.2-B),  employing  only  06  elements 
and  57  node  points,  was  developed  for  interactive  use  on  the 
computer.  This  mesh  behaved  roughly  like  the  finer  meshes 
but  could  be  run  more  frequently,  using  less  computer  time. 

The  goal  in  modelling  the  annulus  was  to  simulate  as 
closely  as  possible  the  actual  arrangement  of  the  lamellae. 


ii mmm 


224  Elements/133  Nodes 
Fig  4.2-A-  Original  Mesh, 


86  Elements/57  Nodes 
Fig  4.2-B.  Coarse  Mesh. 
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The  number  of  fiber  layers,  or  plies,  in  the  annulus  is  esti¬ 
mated  to  be  10  to  12  (Ref  2).  Twelve  plies  were  modelled  in 
the  final  mesh,  bringing  the  total  number  of  elements  to  200 
with  165  node  points  (Fig  4.2-D). 

In  order  to  determine  whether  a  less  refined  mesh  in 
the  annulus  could  be  satisfactorily  used  to  model  its 
behavior,  an  eight-ply  mesh  was  tried  with  252  elements  and 
140  node  points  (Pig  4.2-C).  It.  was  found  that  in  order  to 
match  the  expeci  ment.a  1  initial  displacement  with  the  eight- 
plv  model,  a  slightly  smaller  elastic  modulus  in  the  nucleus 
was  required.  For  example,  for  v  =  .495,  E  was  determined 
to  be  1.62  Kp/sq  cm  (compared  with  1.58  Kp/sq  cm  for  the 
12-ply  model).  Otherwise,  the  two  models  behaved  in  a  very 
similar  fashion. 

In  progressing  from  the  4-ply  model  to  the  12-ply 
model,  the  mesh  size  is  systema fc ica 1 ly  reduced  in  the  annu¬ 
lar  region.  Each  refined  mesh  represents  a  further  sub¬ 
division  within  the  region  and  therefore  a  more  accurate 
solution.  The  similarity  in  behavior  of  the  three  models 
indicates  that  the  solution  is  converging.  Since  the  8-ply 
and  12-ply  models  use  virtually  the  same  amount  of  computer 
time,  it  was  decided  to  use  the  more  physically  accurate 
12-ply  rmde  1 . 

4 . 3  Initia 1  Results :  Visco-Elastic  Nucleus /Elastic  Annulus 

As  noted  in  Section  1.4,  the  disc  was  initially 
modelled  with  a  visco-elastic  nucleus  and  an  elastic 
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12-]>1v  Annulus 


annulus.  Separate  runs  were  made  with  Poisson's  ratio  equal 
to  0.495  and  0.  335.  Msinj  the  trial  and  error  method  pre¬ 
viously  described,  appropriate  values  of  F  wore  found  for 
each  value  of  v  that  allowed  the  initial  displacement  of  the 
model  to  mutch  the  experimentally  observed  displacement. 

To  determine  the  visco-elastic  parameters,  Hinrichsen's 
final  value's  were  chosen  as  a  starting  point; 
qg  =  .01  Kp/sq  cm,  q^  =  27000  Kp-soc/sq  cm.  The  resulting 
curves,  shown  in  Figs  4. 3-A  and  4.3-B,  contrast  sharply  with 
the  expe r inenta 1  curve,  which  rises  smoothly  from  0.014  cm 
to  0.02B  cm  in  20,000  seconds.  For  the  run  in  which  v  was 
set  equal  to  0.333,  the  curve  levels  off  after  500  seconds 
and  remains  constant  at  approximately  0.019  cm.  The  curve 
for  which  *■  was  set  equal  to  0.  495  appears  to  be  oven  more 
eccentric,  rising  only  .0001  cm  in  15,000  seconds.  (The 
roughness  of  this  curve  is  not  significant.  It  is  due  to 
the  expanded  interval  on  the  Y-axis  caused  by  computer 
graphics. ) 

In  order  to  make  these  curves  match  the  experimental 
curve,  the  method  of  varying  qo  and  qj  was  employed. 

Previous  analysis  {Hef  10)  showed  that  the  maximum  displace¬ 
ment  increases  when  qg  is  decreased  and  that  the  slope 
increases  with  decreasing  q^.  Therefore,  the  value  of 
qn  was  progressively  decreased  from  0.01  to  0.0  Kp/sq  cm 
(the  zero  value  corresponding  to  a  Maxwell  Fluid).  The 
vilu“  of  q]  was  varied  from  17000  to  47000  Kp-sec/sq  cm. 

The  results,  [(lotted  in  Fig  4.  3-C  for  v  =  .  3  3  3  ,  show  that  the 
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Displacement  vs  Time,  Elastic  Nucleus 
Variation  of  Visco-elastic  Parameters 


curves  were  unaffected  by  changes  in  qo <  regardless  of  the 
value  of  q^.  Similar  results  were  obtained  with  v  =  .495. 

No  matter  how  the  q  parameters  were  varied,  the  maximum 
displacement  that  could  be  achieved  by  the  jnodel  was 
.0195  an.  In  other  words,  the  overall  disc  model  appeared 
to  he  too  stiff  to  allow  further  displacement. 

Several  approaches  to  the  problem  were  considered.  One 
option  was  to  change  the  values  of  the  orthotropic  material 
constants  in  the  annulus.  However,  these  constants  (with 
the  exception  of  V23  and  v\i )  were  based  upon  previous 
experimental  work.  Even  if  inaccurate,  the  means  were 
lacking  to  determine  seven  additional  constants.  The  value 
of  v?_3  (=  V32)i  on  the  other  hand,  had  been  chosen 
arbitrarily  to  satisfy  Equation  4-3.  Therefore,  a  range  of 
v>23  values  were  tried,  from  0.31  to  0.45.  The  effect  on 
the  (maximum  displacement  was  negligible. 

At  this  point,  an  entirely  different  approach  was 
taken.  Since  the  stiffness  of  the  annulus  appeared  to  be 
preventing  further  deformation  of  the  disc,  it  was  decided 
to  determine  how  "soft"  the  material  constants  in  the  annu¬ 
lus  would  have  to  be  made  to  match  the  experimental 
displacements.  Ry  stipulating  that  the  ratio  of  E3  to  E2  i-n 
the  annulus  remain  constant,  the  value  of  E 3  could  be  varied 
without  affecting  the  values  of  \>]_2’  v13»  etc.  (reference 
Eq  4-1).  Recalling  that  the  nucleus  had  been  observed 
to  behave  hydrostatically  under  axial  load,  it  was  assumed 
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to  hove  the  some  hulk  nodulus  as  water. 


A  Poisson's  ratio 


of  .313  was  chosen  to  simulate  the  deformable,  gel-like 
nucleus  material.  It  was  reasoned  that  if  the  value  of 
Eg  required  to  'natch  the  experimental  displacement  were 
unreasonably  small,  it  could  be  concluded  that  the  nucleus 
as  well  as  the  annulus  must  be  modelled  visco-elastical ly  to 
match  experimental  results.  This  did  indeed  prove  to  be  the 
case.  The  value  of  Eg  was  varied  from  4350  to  .435  Kp/sq 
cm.  The  resulting  maximum  displacements  were  far  below  the 
experimental  value. 

based  upon  these  findings,  it  was  decided  to  model  the 
entire  intervertebral  disc  visco-e 1 ast i cal ly .  The  final 
results  of  this  analysis  are  presented  in  the  following 
chapter . 

4 . 4  Final  Model:  Vise q -E Last ic  Disc 

Once  l he  decision  was  made  to  model  both  the  nucleus 
and  the  annulus  visco-elasticasl ly,  the  remainder  of  the 
investigation  proceeded  smoothly.  Values  of  E  in  the 
nucleus  had  already  been  determined  for  both  v  =  .333  and 
v  -  .495.  It  remained  to  find  the  values  of  <m  and  qg  in 
the  nucleus  and  annulus.  It  should  be  pointed  out  at  this 
point  that  visco-elastic  constitutive  equations  have  yet  to 
be  derived  for  orthotropic  materials.  The  annulus  is  effec¬ 
tively  being  modelled  here  as  orthotropica 1 Ly  elastic  but 
isotropical ly  visco-elastic. 
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Following  the  pioodure  previously  outlined,  the  visco¬ 
elastic  parameters  were  varied  untiL  the  model  matched  the 
experimental  displacement  versus  time  curve.  Fig  4.4-A  shows 
some  intermediate  results  obtained  by  varying  qg  and  qj_  as 


well  as  the  ti 

represent  the 

na)  match.  The  values 

data  are  shown  below. 

of  qo  and  q]_ 

which  best 

V 

no 

<U 

(  nuc  1  eu.s  ) 

(nucleus,  annulus) 

( nucleus , 

annulus  ) 

0.333 

l .  2 

85000 

0.495 

1.5 

85000 

Fig  4 . 4-13 

Note  that  the  same  visco-elastic  constants  were  used  in 
both  the  nucleus  and  annulus.  This  is  most  likely  untrue, 
since  the  materials  in  each  region  are  different 
visco-el ast  ica 1  ly .  However,  it  is  not  possible  to  improve 
upon  these  values  with  the  method  at  hand.  The  problem  is 
that  we  have  four  unknowns  (the  q  values  in  the  annulus  and 
nucleus)  and  only  one  equation;  the  experimentally  observed, 
one-dimensional  displacement  of  the  disc.  Therefore, 
although  we  might  find  other  sets  of  q  values  by  trial  and 
error  which  produce  a  curve  match,  there  would  be  no  way  to 
determine  which  set  is  the  more  accurate.  As  a  matter  of 
curiosity,  other  sets  wore  tried.  Not  surprisingly,  it  was 
found  that  the  experimental  curve  could  be  matched  by  a 
variety  of  q  values.  For  example,  using  the  coarse  mesh  and 
a  v  value  of  .333,  the  curve  was  matched  by  the  following 
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Visco-Elastic  Nucleus  and  Annulus. 


pa  rameters : 


Nucleus  Annulus 


do 

dl 

d0 

dl 

Bet 

n 

0.8 

81000 

0.8 

81000 

Set 

#2 

1.0 

70000 

0.  6 

90000 

In  other  words,  the  experimentally  observed,  one¬ 
dimensional  behavior  of  the  vertebra  could  be  accurately 
simulated  by  the  model  using  widely  varying  values  for  qo / 
q i ,  and  v.  The  q  values  which  were  chosen  to  represent  the 
final  model  are  those  given  in  Fig  4.4-B. 

4 . 5  Displacement  Prof i les 

In  order  to  visualize  the  deformations  taking  place 
throughout  the  joint,  as  time  increases,  plots  of  the 
deformed  model  were  made  at  t  =  0,  t  =  70,  and  t  =  170 
minutes.  Separate  plots  were  generated  using  Poisson's 
ratios  in  the  nucleus  of  v  =  .333  and  v  =  .495.  Since  the 
resulting  plots  were  virtually  identical,  we  conclude  that 
the  displac  ements  do  not  depend  on  the  value  of  v  in  the 
nucleus  as  long  as  it  is  paired  with  an  appropriate  E  value. 
For  this  reason,  only  the  v  =  .333  plots  are  presented  here. 

Figure  4.5-A  presents  the  undeformed  shape  of  the  model. 
Figures  4. 5-B,  C,  and  D  are  the  deformed  shapes  at  t  =  0, 
t  =  70,  and  t  =  170  minutes,  respectively.  To  more  clearly 
depict,  the  overall  changes  taking  place,  Fig  4.5-E  shows  the 
t  =  170  plot,  overlaid  on  the  undeformed  (t  =  0-)  plot,  with 
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Fin  -1.5-C.  Deformed  Diet,  t  =  70  Minutes. 


Fiq  4.5-D. 


Deformed  Plot,  t 


170  Minutes. 


Fict  4.5-E.  Comparison  of  Unde  formed  Shape  to 
Deformed  Shape  at  T  =  170  Minutes. 
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the  mesh  lines  deleted.  Finally,  for  the  purpose  of 
comparison,  a  deformed  plot  of  Hinrichsen's  homogeneous 
model  (Ref  10)  at  t  =  1 R 0  minutes  is  reproduced  in  Fig 
4.5-F.  In  each  of  these  figures,  the  deformations  have  been 
magnified  so  that  they  may  be  more  readily  observed.  The 
magnification  factors  are  noted  in  the  figures. 

Note  that  there  are  striking  differences  between  the 
behavior  of  Flinrichsen  ’  s  homogeneous,  isotropic  disc  and 
that  of  the  present  model.  The  homogeneous  disc  spreads 
outwardly  in  the'  radial  direction  with  time.  Its  radius  at 
the  centerline  increases  by  approximately  4%  over  180  min¬ 
utes  of  loading.  On  the  other  hand,  the  orthotropical ly 
modelled  annulus  exerts  a  considerable  restraining  force  on 
the  outward  flow  of  the  nucleus.  After  a  similar  period  of 
time,  the  radius  of  the  disc  at  the  centerline  increased 
only  1%. 

Within  the  non -homogeneous  disc,  the  annulus  can  be 
seen  to  be  adjusting  to  the  outward  pressure  of  the  nucleus. 
.Since  the  annular  fibers  are  firmly  attached  to  the  end 
plates,  they  are  constrained  from  radial  movement  at  the 
attach  points.  Therefore,  the  radial  deformation  of  the 
annulus  is  greatest  at  the  centerline  of  the  disc.  Note 
that  the  fiber  layers  nearest  to  the  annulus  are  more 
directly  affected  by  its  outward  flow.  The  inner  ply  con¬ 
sequently  displace  farther  radially  than  do  the  outer  ply. 
Hence,  one  would  expect  to  see  higher  compressive  stresses 
in  the  inner  ply  than  in  the  outer  ply.  The  attachment  of 
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the  annular  fibers  to  the  end  plate  should  have  an  effect  on 
the  stresses  within  the  endplate  as  well,  particularly  on 
the  shear  stress. 

The  vertebral  body  and  the  bony  end  plate  appear  to 
undergo  virtually  no  deformation  with  time.  The  entire 
axial  displacement  takes  place  within  the  intervertebral 
disc.  This  is  not  surprising  because  of  the  large  dif¬ 
ference  between  the  elastic  modulus  of  bone  versus  that  of 
f  ibroca r  t i  lage . 

In  the  following  section,  the  actual  stress  redistribu¬ 
tions  are  presented  and  compared  to  those  produced  in  the 
homogeneous  disc  model. 

4  .  f>  Stress  Redistributions 

Horizontal  profiles  of  the  radical,  axial,  hoop,  and 
shear  stress  components  were  taken  at  three  different  joint 
levels;  through  the  disc,  the  bony  end  plate,  and  the  ver¬ 
tebral  body.  Fig  4.5-A  indicates  the  position  of  these 
horizontal  slices.  To  determine  the  effect  of  the  Poisson's 
ratio  in  the  nucleus,  separate  profiles  were  generated  for 
v  =  .333  and  v  =  .495.  It  was  found  that  for  the  most  part, 
the  v  value  had  little  effect  upon  the  behavior  of  the 
model.  A  few  differences  were  noted  in  the  way  that  certain 
stress  components  varied  with  time  in  the  disc  and  end 
plate.  For  example,  with  v  =  .333,  the  hoop  and  radial 
stresses  in  the  nucleus  increase  approximately  70%  over  a 
170  minute  interval.  The  same  stresses  decrease  10%  during 
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that  interval  wnen  v  -  .495.  These  differences  may  be 
attributed  to  the  fact  that  the  higher  value  of  v  decreases 
the  rate  of  displacement  in  the  nucleus  by  limiting  the 
direct  strain.  At  this  point  a  choice  had  to  be  made  as  to 
which  value  of  v  would  be  used  in  the  final  model.  Either 
value  could  be  used  to  match  the  experimental  displacements 
and  both  produced  similar  stress  profiles.  It  was  decided 
that  the  lower  value,  corresponding  to  a  deformable,  gel¬ 
like  solid,  would  provide  the  best  model  of  the  behavior  of 
the  nucleus.  The  resulting  horizontal  stress  profiles  are 
presented  in  Figs  4.6-C  through  4.6-N  for  times  varying  from 
0  to  170  minutes. 

A  second  set  of  profiles  was  taken  to  show  the  stress 
redistributions  in  the  individual  lamella.  Three  separate 
plies  are  profiled;  an  inner  ply ,  an  outer  ply,  and  one  mid¬ 
way  between  (see  Fig  4.6-B).  The  stresses  are  averaged  over 
the  length  of  the  ply  and  plotted  against  time.  To  compen¬ 
sate  for  the  effect  of  the  alternating  fiber  angle,  the 
stresses  ir.  two  adjacent  lamella  were  initially  averaged 
together  to  represent  the  stress  in  a  single  ply.  As  it 
turned  out,  the  stresr.es  vary  so  smoothly  from  ply  to  ply 
that  this  averaging  technique  proved  unnecessary. 

Therefore,  the  stress  profiles  presented  in  Figs  4.6-0 
through  4.6-R  represent  stresses  averaged  within  single 


plies. 


We  first  examine-  the  stress  redistributions  in  the 
disc,  where  the  effect  of  the  annu ins-nucleus  interface  is 
immediately  evident.  At  t=0  minutes,  both  the  hoop  and 
axial  stresses  decrease  71%  across  the  interface 
(Figs  4.6 -n,  E).  The  shear  stress,  which  is  virtually  zero 
in  the  nucleus,  goes  to  .2  Kp/sq  cm  in  the  annulus 
(Fig  4.6-F).  As  time  progresses,  and  the  disc  adjusts  to 
the  load,  these  changes  become  loss  abrupt.  Comparing  these 
profiles  to  similar  ones  made  by  Hinrichson  (Ref  10)  of  a 
homogeneous  disc,  we  see  that  the  stress  magnitude  and  the 
manner  in  which  they  vary  in  time  are  very  similar.  The 
general  shapes  of  the  profiles  arc  also  the  same,  although 
the  effect  of  the  annulus  is  necessarily  absent  in  the  homo¬ 
geneous  model. 

Looking  at  the  physical  structure  of  the  annulus,  one 
would  expect  its  behavior  to  be  that  of  a  thick-walled 
cylindrical  shell  subjected  to  internal  pressure.  That  is, 
one  expects  to  see  dominant  tensile  hoop  stresses  in  the 
annulus  as  a  result  of  axial  loading  and  outward  pressure  of 
the  nucleus.  Kulak  (Ref  15)  demonstrated  that  this  is,  in 
fact,  the  case  in  lumbar  discs.  On  the  other  hand,  thoracic 
discs,  such  as  the  specimen  modelled  in  this  analysis,  are 
relatively  thinner  and  therefore  more  constrained  by  their 
end  plates.  Figs  4.6-E  and  4.6-Q  show  that  the  annulus  is 
in  a  state  of  hoop  compression.  It  is  acting  in  effect,  not 


like  a  shell,  hut  like  a  ring  umlcr  axial  load.  The  magni¬ 
tude  of  this  compression  is  small;  approach  ing  zero  at  the 
outer  edge  of  the  annulus.  This  agrees  with  Kulak's  analy¬ 
sis  of  thoracic  discs.  He,  too,  observed  small  hoop 
stresses  in  the  annulus  and  compression  in  many  of  the  annu¬ 
lus  fibers. 

Fig  4.6-C  shows  that  the  radial  stress  tends  to  vary 
smoothly  across  the  interface  at  all  times,  though  it 
decreases  rapidly  in  th>>  outer  ply  of  the  annulus.  This 
indicates  that  the  outer  plies  are  undo),  considerably  less 
radial  compress  ion  than  the  inner  plies.  The  same  conclu¬ 
sion  may  be  drawn  from  Fig  4. 6-0  which  shows  36%  more  radial 
compression  in  the  inner  ply  than  in  the  outer.  Note  that 
in  the  outer  ply,  all  four  stress  components  are  fairly 
constant  over  time  (Figs  4. 6-0, P,Q, R) .  In  contrast,  the 
stresses  in  the  inner  ply  (with  the  exception  of  shear)  tend 
to  increase’  with  time.  This  implies  that  the  inner  lamellae 
do  most  of  the  work  in  restraining  the  nucleus  under  small 
loads . 

In  the  bony  end  plate,  the  magnitude  of  the  shear, 
hoop,  and  radial  stress  components  is  much  greater  than  in 
any  other  region  of  the  joint.  This  signifies  the  impor¬ 
tance  of  the  end  plate  as  a  structural  member  of  the  joint. 
Although  the  displacement  profiles  indicate  that  little 
tctual  deformation  takes  place  in  the  end  plate  under  axial 
load,  it  do<s  undergo  significant  redistributions  of  stress 
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over  time.  For  example,  Fig  4.6-1  shows  that  the  hoop 
stress  in  the  end  plate  increases  16%  in  a  170  minute 
interval.  In  contrast,  the  radial  stress  increases  with 
time  only  in  the  region  directly  above  the  nucleus  (Fiq 
4.6-G).  The  axial  an  1  shear  stresses  in  the  end  plate  also 
appear  to  l>e  affected  by  their  proximity  to  the  annulus- 
nucleus  interlace  (Figs  4.6-11,  4.6-d).  It  is  unlikely  that 
the  abrupt  peak  in  the  axial  stress  curve  at  r  ■=  i.  9  cm  is  a 
completely  accurate  representa t ion  of  the  stress  in  that 
area.  The  peak  may  be  partially  due  to  insufficient  refine¬ 
ment  of  the  mesh  size  in  that  region.  However,  further 
refinement  of  the  mesh  would  bo  very  expensive  in  terms  of 
computer  resources.  Also,  it  is  questionable  whether  any 
degree  of  ref iixronl  would  produce  a  smooth  stress  tran¬ 
sition  in  a  region  whore  three  different  types  of  materials 
are  joined. 

In  the  vertebral  body,  the  hoop  and  shear  stresses 
(Fiqs  4.6-M  and  4.6-N)  are  virtually  constant  across  the 
radius  and  do  not  vary  with  time.  Note  that  the  last  data 
point  on  these  plots  represents  the  stress  in  the  thin  shelL 
of  cortical  hone  which  encases  the  vertebral  core.  The 
stresses  in  the  shell,  like  those  in  the  bony  end  plate,  did 
vary  with  time.  It  is  significant  to  note  that  the  cortical 
shell  is  in  a  state  of  hoop  tension.  That  is,  its  boh.v’or 
is  that  of  a  classical  shell.  In  both  fho  vertebral  l>>dy 
and  the  end  plate,  the*  axial  stress  tends  to  decrease  with 
radius,  but  dramatically  increases  at  the  outer  edge  of  the 


joint.  This  supports  the  idea  that  the  load  is  transmitted 
mainly  through  the  cortical  shell  and  end  plate  to  the 
nucleus,  which  then  acts  as  the  primary  load  carrying  struc¬ 
ture  of  the  joint. 
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V. 


CONCLUSION 


The  mechanical  behavior  of  the  intervertebral  joint  can 
be  realistically  modelled  through  the  use  of  the  finite  ele¬ 
ment  method.  By  matching  one-dimensional  experimental 
displacement  data  material  parameters  are  determined  which 
provide  a  good  overall  description  of  the  time  dependent 
behavior  of  the  joint. 

The  inclusion  of  the  orthotropic  material  properties  of 
the  annulus  in  the  model  caused  significant  stress  redistri¬ 
butions  throughout  the  disc.  These  redistributions  indicate 
the  importance  of  the  annulus  as  a  restraining  mechanism  on 
the  outward  flow  of  the  nucleus  under  axial  loading.  It  is 
concluded  that  the  orthotropic  'nature  of  the  annulus  is  a 
factor  which  should  be  included  when  modelling  the  interver¬ 
tebral  joint. 

The  results  of  this  study  indicate  that  both  the 
nucleus  pulposis  and  the  annulus  fibrosis  are  visco-elastic 
materials.  The  initial  assumption  that  the  experimentally 
observed  creep  response  of  the  joint  was  due  to  the  nucleus 
alone  proved  unworkable.  It  was  found  that  the  creep  beha¬ 
vior  of  the  joint  could  be  simulated  only  if  both  the 
nucleus  and  the  annulus  were  modelled  visco-e last ica 1 ly . 

Due  to  the  limitations  of  the  curve  matching  technique, 
unique  values  for  the  visco-elastic  parameters  in  these 
regions  cannot  be  determined.  Determination  of  unique 
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parameters  would  require  at  least  two-dimensional  displace¬ 
ment  data,  which  at  present  is  unavailable. 

Examination  of  the  displacement  profiles  and  stress 
redistributions  within  the  model  provides  insight  into  the 
internal  behavior  of  the  joint  under  axial  loading.  It  was 
found  that  the  cortical  bone  which  encases  the  vertebral 
body  behaves  essentially  like  a  shell  subjected  to  internal 
pressure.  On  the  other  hand,  the  geometry  of  the  thoracic 
disc  is  such  that  the  annulus  behaves  like  an  axially  loaded 
ring.  Under  small  Loads,  the  inner  lamellae  of  the  annulus 
do  most  of  the  work  in  restraining  the  outward  flow  of  the 
nucleus.  It  appears  that  when  the  joint  is  subjected  to 
axial  compression,  the  load  is  transmitted  mainly  through 
the  cortical  shell  and  end  plates  to  the  nucleus,  which  then 
acts  as  the  primary  load  carrying  structure  of  the  joint. 
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^ An  understanding  of  the  mechanical  properties  and  behavior  of  the 
intervertebral  disc  is  critical  to  several  current  areas  of  research. 

Among  these  are  the  study  of  the  effects  of  extreme  gravitational  forces  on 
the  air  crews  of  high  performance  aircraft  and  the  related  problem  of  spi¬ 
nal  injuries  due  to  aircraft  ejection.  This  study  was  undertaken  in  order 
to  construct  a  realistic  analytical  model  of  the  intervertebral  joint  using  - — - 

the  finite  element  approach. 
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Rxper iments  Iwvo  -.howi  that  healthy  Intervertebral  joints  exhibit 
creep  when  subjected  to  axial  load.  The  focus  of  this  investigation  is  the 
behavior  and  interaction  of  the  nucleus  and  annulus  under  loading.  An  axi- 
symmetrlc  finite  element  model  is  employed  which  incorporates  a  linear 
visco-elagtie  constitutive  relation  for  the  annulus  and  nucleus,  based  upon 
a  three-parameter  Kelvin  solid.  The  visco-eiastic  constants  are  found  by 
matching  one-dimensional  axial  experimental  data  with  this  two  dimensional 
model.  The  nucleus,  which  is  composed  of  a  series  of  concentric  lamellae 
of  collagen  fibers,  has  been  modelled  as  a  12  nlv  structure  of  orthotropic 
material . 

Results  are  presented  which  depict  the  displacement  profiles  and 
stress  redistributions  occurring  as  a  consequence  of  the  interaction  of  the 
annulus,  nucleus,  and  bony  end  plate  under  axial  compressive  loading. 

These  profiles  clearly  show  the  importance  of  the  nucleus  as  a  load 
carrying  structure  uni  the  action  of  the  annulus  in  restraining  the  outward 
flow  of  the  nucleus.  Also  presented  are  the  variations  in  stress  with  time 
in  three  of  the  lamina  o‘  the  annulus.  Results  indicate  that  the  orthotro¬ 
pic  material  properties  of  the  annulus  have  a  significant  impact  upon  its 
time  dependent  behavior  an  1  should  be  included  when  modelling  the  joint. 


